Dynamical mean-field theory for light fermion-heavy boson mixtures on optical lattices 
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We theoretically analyze Fermi-Bose mixtures consisting of light fermions and heavy bosons that are loaded 
into optical lattices (ignoring the trapping potential). To describe such mixtures, we consider the Fermi-Bose 
version of the Falicov-Kimball model on a periodic lattice. This model can be exactly mapped onto the spinless 
Fermi-Fermi Falicov-Kimball model at zero temperature for all parameter space as long as the mixture is ther- 
modynamically stable. We employ dynamical mean-field theory to investigate the evolution of the Fermi-Bose 
Falicov-Kimball model at higher temperatures. We calculate spectral moment sum rules for the retarded Green's 
function and self-energy, and use them to benchmark the accuracy of our numerical calculations, as well as to re- 
duce the computational cost by exactly including the tails of infinite summations or products. We show how the 
occupancy of the bosons, single-particle many-body density of states for the fermions, momentum distribution, 
and the average kinetic energy evolve with temperature. We end by briefly discussing how to experimentally 
realize the Fermi-Bose Falicov-Kimball model in ultracold atomic systems. 

PACS numbers: 03.75.Lm, 37.10.Jk, 67.85.-d, 67.85.Pq 



I. INTRODUCTION 

Experimental work in ultracold atomic systems in optical 
lattices has been progressing rapidly. First generation exper- 
iments focused primarily on single species systems, or mix- 
tures of different isotopes of the same atomic species. Now, 
experimental work is progressing into the realm of mixtures 
of different species of atoms. Since these atomic species of- 
ten have significantly different masses, one immediately ex- 
pects such systems to respond differently than isotopic mix- 
tures. In addition, if light alkali atoms like Li or K are used, 
or if one uses alkaline earth atoms like Sr or rare-earths like 
Yb, then one has the prospect for modifying the particle statis- 
tics from Fermi-Dirac to Bose by simply changing the isotope 
employed in the experiment. 

In this work, we focus on Fermi-Bose mixtures where the 
bosonic atoms are the heavy atoms (which is typically the 
experimental situation if Li or K is the light fermion and 
Rb, Cs, Sr, or Yb is the heavy boson). While much work 
on Fermi-Bose mixtures has focused on either how the pres- 
ence of the fermions modifies the Bose-Einstein condensa- 
tion (BEC) of the bosons, or how the presence of the bosons 
modifies the interactions of the fermions {e.g. by allowing 
phonons) QJJ] , our work here focuses on the opposite situa- 
tion where the heavy bosons are so heavy that we can ignore 
the quantum-mechanical effects of their kinetic energy, and 
hence they never condense into a BEC. Instead, at low tem- 
perature, the mixture either phase separates or forms static 
density wave patterns, which can be quite complex. We be- 
lieve such a situation arises when the tunneling amplitude of 
the heavy bosons on the optical lattice is more than an or- 
der of magnitude smaller than the tunneling amplitude of the 
light fermions. This typically occurs when the lattice depth 
is deep (more than 10 to 15 recoil energies of the Rb for a 
K-Rb mixture) and the system is well represented by a single- 
band model yfl- In this situation, experiments are likely to 
be run at temperatures significantly higher than the BEC tem- 
perature, and it is more likely that static density waves would 



form instead of superfluidity (in other words, we are exam- 
ining how the fermions modify the Mott-insulating state of 
the bosons, which occurs when the boson-boson repulsion is 
much larger than the boson tunneling amplitude [3]). Long- 
range effective boson-boson interactions are generated via the 
interactions with the mobile fermions. This realm is not so 
well known within the atomic physics community, although 
similar models have been widely studied within the condensed 
matter physics community, as we describe below. 

Recently, mixtures of fermionic 40 K and bosonic 87 Rb 
atoms have been first studied at fixed interspecies interaction 
strengths by two different experimental groups Q,§], and later 
with tunable interactions |6!], where a shift of the bosonic su- 
perfluid to Mott insulator transition has been observed due to 
the interspecies coupling (irrespective of its sign but with sig- 
nificant asymmetry). There are several theoretical proposals 
to explain this effect 0, H, H, OH El] . Motivated by these ex- 
periments, here we study such light-Fermi-heavy-Bose mix- 
tures with the Fermi-Bose version of the Falicov-Kimball 
(FK) model |fi2t[T3l[Tl, the Fermi-Fermi version of which 
has been widely discussed in the condensed matter litera- 
ture d M- 

In our case, the bosons have no quantum dynam- 
ics of motion, but they can sample all possible heavy-atom 
configurations in an annealed statistical-mechanical sense. 

Our main results are as follows. First, we examine the sym- 
metries of the Hamiltonian, and show that the Fermi-Bose FK 
model can be mapped exactly onto the spinless Fermi-Fermi 
FK model at zero temperature for all parameter space as long 
as the mixture is thermodynamically stable. Since this map- 
ping is only approximate at low temperatures and it fails at 
high temperatures, we develop dynamical mean-field theory 
(DMFT) IU7I1 (which becomes exact in infinite dimensions) to 
investigate the effects of temperature and how the Fermi-Bose 
system evolves into an effective Fermi-Fermi system. In addi- 
tion, we calculate spectral moment sum rules for the retarded 
Green's function and self-energy, and use them to check the 
accuracy of our numerical calculations, as well as to reduce 
the computational cost IU8l fl9L uOl 12111 . We also present typ- 
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ical numerical results for the Fermi-Bose FK model includ- 
ing the occupancy of the bosons, single-particle many-body 
density of states (DOS) for the fermions, momentum distri- 
bution, and the average kinetic energy. The DOS has signifi- 
cant temperature dependence as the system evolves from high 
to low temperatures. For example, in the insulating regime, 
where the magnitude of the boson-fermion interaction is much 
larger than the magnitude of the fermion tunneling amplitude, 
many 'upper Hubbard bands' exist at high temperature, with 
the weights of each band determined by the probability that 
the bosons multiply occupy a lattice site. These bands then 
evolve to just one upper and one lower Hubbard band at low 
temperature. 

The remainder of this manuscript is organized as follows. 
In Sec. UH we first introduce the Hamiltonian for the Fermi- 
Bose version of the FK model, and then discuss its symme- 
tries. In Sec. [HI] we develop the DMFT for this model, where 
we also discuss spectral moment sum rules for the retarded 
Green's function and self-energy. The numerical results ob- 
tained from the DMFT are discussed in Sec. [IV] and a brief 
summary and conclusions are presented in Sec.lVl 



II. LIGHT FERMION-HEAVY BOSON MIXTURES 

We analyze Fermi-Bose mixtures consisting of light 
fermions and heavy bosons that are loaded into optical lat- 
tices. We consider the limiting case such that the band mass of 
bosons Mf, is much greater than that of fermions (M& 3> Mf). 
In this case, the tunneling term for the bosons tb is much 
smaller than that of the fermions (tb <C tf), and we may set 
tb = in the quantum-mechanical Hamiltonian that describes 
the system. Such mixtures could be experimentally realized in 
ultracold atomic systems. For instance, a 40 K- 87 Rb mixture 
with M b m 2.2M/, a 6 Li- 41 K mixture with M b « 6.8M/, 
or a 6 Li- 133 Cs mixture with Mb ~ 22. 2Mj are good candi- 
dates for the applicability of this Hamiltonian. In addition, one 
could also create species-dependent optical lattices for differ- 
ent isotopes of the same atom such that the bosonic isotope is 
localized but the fermionic one is not, which might be most 
relevant for isotopic mixtures of Sr or Yb (in other words, 
it isn't the masses per se that need to be different, but it is 
the tunneling amplitudes when in an optical lattice that need 
to be quite different). It is important to emphasize that, al- 
though we consider the limit of localized bosons, the problem 
at hand is still a strongly correlated many -body problem since 
the bosons and fermions are coupled by the interaction and we 
take annealed statistical averages over all heavy-atom config- 
urations. 



A. The Fermi-Bose Falicov-Kimball (FK) Model 

With such systems in mind, we analyze the following 
Hamiltonian to describe the light fermion-heavy boson mix- 
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where /• (fi) is the creation (annihilation) operator for an itin- 
erant spinless fermion at site i, and b\ (bi) is the corresponding 
operator for a localized boson at site i. The fermionic op- 
erators satisfy the usual canonical anticommutation relations 
{/i> f}} — $ij where 5ij is the Rronecker delta function, 
the bosonic operators satisfy the usual canonical commuta- 
tion relations = <%, and the fermionic and bosonic 
operators all mutually commute. The first term is the ki- 
netic energy of the fermions with tf denoting the nearest- 
neighbor tunneling. The second and third terms are the on-site 
density-density interactions between the bosons themselves, 
and between the bosons and fermions, respectively. The on- 
site fermion-fermion interaction is not allowed because of the 
Pauli exclusion principle, i.e. Uff — > oo. The last two 
terms involve the chemical potentials of the fermions (/i /) and 
bosons (fib), that are employed to adjust the filling of the cor- 
responding particles to the desired values. 

The model Hamiltonian given in Eq. (Q]) then corresponds to 
the Fermi-Bose version of the FK model, the Fermi-Fermi ver- 
sion of which has been studied extensively in the condensed- 
matter literature lfl5l [T^l . At low temperatures, the Fermi- 
Fermi version of the model is known to possess significant 
regions of density wave order with complicated patterns lfl5ll . 
and also has a strong tendency to phase separate when the ligh t 
fermion-heavy fermion interaction is large and repulsive II22I1 . 
It is also known that the heavy fermions cannot generate an 
effective retarded light fermion-light fermion attraction lead- 
ing to superconductivity as a consequence of Anderson's the- 
orem 112311 - Within the DMFT context, essentially all response 
functions (both static and dynamical) and all kinds of trans- 
port have been evaluated for the system [16]. We show below 
that many of these low temperature properties are shared by 
the Fermi-Bose FK model, but at higher temperatures, the be- 
havior is quite different. 

B. Symmetries of the Hamiltonian 

This Hamiltonian possesses partial particle-hole symme- 
tries for the fermions and bosons. These symmetries hold on 
a bibartite lattice where it is possible to divide the entire lat- 
tice into two sublattices A and B such that the fermionic tun- 
neling only connects different sublattices. When the particle- 
hole transformation is applied to the fermions, i.e. fj — > 
// 1 (-1)p« and ft -> {ff)H-l) p(l) withp(i) = fori e A 
and p(i) = 1 for i G B, then, up to a numerical shift, it can be 
shown that the Hamiltonian that is expressed in terms of the 
hole operators for the fermions maps onto the starting Hamil- 
tonian with fif — ► — [if, Ubf — > — Ubf and fib —* Hb — Ubf- In 
the canonical ensemble, this means that the energies are sim- 
ply related by E(p f ,p b , U bf ) = E(l - pf,p b , -Ubf), where 
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Pf and pb are the fillings for fermions and bosons, respec- 
tively. We notice that this mapping for the fermions holds at 
any temperature T, unlike that for the bosons, discussed next. 

When the particle-hole transformation is applied to the 
bosons, i.e. b\ -► b£(-l) p V) and b, -> (&£)+(- 
then, up to a numerical shift, it can also be shown that the 
Hamiltonian that is expressed in terms of the hole opera- 
tors for the bosons maps onto the starting Hamiltonian with 
Pf — > pf + Ubf and pi, — ► pb + Ubb- In the canonical ensem- 
ble, this means that the energies of many-body eigenstates are 
simply related by E(pf,p b , Ubf) = E(pf, 1 + pb, Ubf). This 
is not enough to map Green's functions onto each other. The 
subtle issue is that the particle states with rib = 0, 1, 2, • • • , oo 
map onto the hole states n' b l = rib + 1 = 1, 2, 3, • • • , oo, hence 
we lose all information about the weight of the n' b l = state 
when more than one rib > state has a nonzero density; when 
only two states are nonzero, the mapping above plus the rela- 
tion between particle and hole densities allows us to find all of 
the bosonic weights. This occurs exactly only at T = where 
at most two states have nonzero weights. Hence the bosonic 
particle-hole symmetry only occurs at T = 0. This symmetry 
also holds in the atomic (£& = 0) limit of the Bose-Hubbard 
model at T = 0, where the energy for the particle excitations 
of the ri(,th state is degenerate with that for the hole excita- 
tions of the (rib + l)th state, but it does not hold when tb ^ 0, 
because one can no longer label the many-body eigenstates by 
the local boson particle number. 

The latter symmetry implies that there can be at most two 
bosonic states that can be occupied at T = 0. For instance, 
rib = {0, 1} states for < pb < T, rib = 1 state for pb = 1; 
rib = {1,2} states for 1 < pb < 2; rib = 2 state for pb = 2; 
rib = {2, 3} states for 2 < pb < 3; etc. This is similar to what 
happens in the spinless Fermi-Fermi FK model where there 
can be at most two fermionic states for the heavy fermions due 
to the Pauli exclusion principle. Therefore, the Fermi-Bose 
FK model can be exactly mapped onto the well-studied spin- 
less Fermi-Fermi FK model for all parameter space at T = 
as long as the mixture is thermodynamically stable. Since this 
mapping is only approximate at low temperatures, and it fails 
at high-temperatures, we develop DMFT to investigate the ef- 
fects of temperature. 



HI. DYNAMICAL MEAN-FIELD THEORY 

The DMFT for the Fermi-Fermi FK model is well estab- 
lished in the literature lfl6l [nil , and it can be easily general- 
ized to the Fermi-Bose case that we consider here. For this 
purpose, we notice that the bosons cannot undergo BEC since 
there are a fixed number of localized bosons on each lattice 
site and the local boson particle number is an operator that 
commutes with the Hamiltonian. Then the impurity partition 
function can be integrated analytically since the effective ac- 
tion is quadratic in the fermionic operators and depends only 
on the local bosonic number operators. As a result, similar to 
the Fermi-Fermi case, there are three main equations that need 
to be solved self-consistently for the Fermi-Bose FK model. 



A. Imaginary-axis Formalism 

The first equation is the single-particle lattice Green's func- 
tion for the fermions defined by the time-ordered product 
G(r) = -(T r / 4 (r)/ i t (0)) in the imaginary time (0 < r < f3) 
representation, where T T is the imaginary time ordering op- 
erator, and (O) = Tr{e~" H 0}/Z is the ensemble average 
with Z = Tr{e _/3ff } the partition function and (3 = 1/T 
the inverse temperature. Here, the creation and annihila- 
tion operators are in the imaginary-time Heisenberg repre- 
sentation 0(t) = e HT O(0)e~ HT . In the Matsubara fre- 
quency representation, G(iuj n ) — dTe lu)nT G{t), where 
uj n = (2n + 1)tt/P is the fermionic Matsubara frequency, 
the Green's function becomes 



T^o G^iiuin) - Ubf rib' 



(2) 



where w 7lb = Z n jZ (with Z = £)~ =0 Z nb ) is the probabil- 
ity for each site to be occupied exactly by n& = 0, 1, 2, • • • , oo 
bosons, and Gq 1 (iuJ n ) = rio n + pf — X(iLu n ) is the bare 
Green's function with \{iuj n ) the dynamical mean field. Here, 

Z nb = e Pl^ b -U bb n b (n b -l)/2] Zo{ ^ _ jj^ (3) 

is the partition function for the rib state, where 



Zo(pf)=2e^/* l ^ + ^- X ^ 



lUJ n 



(4) 



is the fermionic partition function for Ubf = 0. The prefactor 
in Eq. is added to give the correct noninteracting result 
when X(ioj n ) = 0. In the case of the spinless Fermi-Fermi 
FK model, due to the Pauli exclusion principle, i.e. the heavy 
fermion-heavy fermion interaction Ubb — ► oo, only two states 
(too and wi) can be occupied such that wi = 1 — wo is the 
filling of localized fermions. Note these equations are similar 
to the solution of the classical Holstein model M24I1 . but here 
the boson states are discrete, while there they are continuous. 

The partition functions also satisfy the well-known rela- 
tions Zo(pf) = DetG^" 1 (ia; n ) and Z = DetG -1 (iui n ), so 
that the bare Green's function can be re-expressed as 



G 1 (iuj n ) = G 1 (iuj n ) + Y, 1 (iuJ n ), 



(5) 



which is our second equation. This is Dyson's equation which 
relates the bare Green's function to the self-energy £(w n ), 
and can also be thought as the definition of the self-energy. 

The third equation is given by the lattice Hilbert transform 
of the noninteracting DOS, 



G(iu> r , 



dep(e) 



(6) 



where G(iu> n ) = ^ k G(k, iu> n ) with G(k, iui n ) — l/[iu n + 
Pf — Yj(ioj n ) — e(k)j the momentum-resolved Green's func- 
tion, and p(e) = £ k <S[e - e(k)] = e^/**) 2 /{^t*) 
is the noninteracting DOS for the infinite-dimensional hy- 
percubic lattice with S(x) the delta function and e(k) = 
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— 2i/ 53<=i cos(fcja) the energy dispersion for the fermions. 
Here, we used the fact that the tunneling tf scales with di- 
mension d such that tf = t*/y/4d J25|. This equation can 
be rewritten in terms of the Faddeeva function G(iu) n ) — 
-i^e-^+Vf-^^ 'erfc{-ir)[iu n + p f - S(io;„)]}, 
where erfc(z) = (2/y/n) J z °° dze~ z is the complex comple- 
mentary error function, and 77 = Sign{Im[icj n — T,(ioj n )]}. 

It is convenient to solve Eqs. ([2), ( a nd © self- 
consistently by using an iterative approach fcdl . For a fixed 
set of C/bf,, Ubf and T, our strategy is as follows: (1) choose 
initial values for pf and p^, (2) start with an initial value for 
Ti(iuj n ), e.g. E(z<j„) = 0, and plug it into Eq. (0 to solve for 
G(iuj n )\ (3) use G(iuj n ) and the initial T,(iui n ) in Eq. (0 to 
obtain Go(iui n )', (4) plug Go(iuj n ) into Eq. (0 to solve for the 
new G(iu) n ); (5) use Go(ioj n ) and the new G(iu> n ) in Eq. (0 
to obtain a new E(iw n ); and (6) iterate steps (2)-(5) with the 
new E(iu/ n ) until the solution converges. Finally, we (7) ad- 
just /if and fij, until the desired pf = (fjfi) and pb = {b\bi) 
values are reached. Typically, the self-energy converges to 
eight decimal points in less than twenty iterations (in a few 
seconds) for each choice of pf and /if,, and the entire proce- 
dure takes less than a minute to obtain the particular values 
for pf and pf,. 



B. Real-axis Formalism 

Once the chemical potentials pf and pb, and the occu- 
pation probabilities w nb are calculated from the imaginary- 
axis calculation, we can employ the analytic continuation of 
Eqs. (0, (0 and (0 with iui n — > u> + i0 + , to calculate G(w) 
and E(w) on the real axis using a similar iterative approach. 
Typically, the convergence is slower on the real axis than on 
the imaginary one, especially near the correlation induced 
band edges which leads to a number of poles in the self-energy 
when Ubf is large enough. Therefore, it is important to verify 
the consistency between the imaginary- and real-axis calcula- 
tions. 

For instance, one of the stringent tests is the compari- 
son of the Green's function that is calculated directly on the 
imaginary axis by the algorithm described above to the spec- 
tral representation for the imaginary axis Green's function 
given by G(iuj n ) = -(l/7r)Im J*^ dwG(uj)/(iu> n - ui), 
with the Green's function generated by the real-axis code 
appearing in the integrand. We usually achieve more than 
four digits of accuracy between the two calculations. An- 
other test is the comparison of the filling of the fermions 
that can be calculated directly on the imaginary axis via 
Pf = {1/(3) Y^=~oo G{iuo n ) to that found from the real axis 
with pf = — (l/V) duF(cj)Im.G(u), where F(u) — 
1/(1 + e' 3 ") is the Fermi-Dirac distribution function. 

One can also check the spectral moment sum rules for 
the retarded Green's function, and also for the retarded self- 
energy 111 8[ Il9l |2(A Kill . These moments are integrals of 
powers of frequency multiplied by the corresponding spec- 
tral function and integrated over all frequency. For instance, 
the spectral moments for the Green's function are defined as 



/j^j = — (l/V) J_ duJUJ m lmG(u)). It can also be shown that 

M « = ({[fi, H] m , //}), where {O u 2 } = O x 2 + 2 1 is 
the anti-commutator, and [Oi,02]m is the multiple commu- 
tation operator such that [0i,O2]i = [Oi,02] = O1O2 — 
O2O1; [Oi,0 2 ] 2 = [[Oi,0 2 ],0 2 ]; etc. Evaluating the com- 
mutators is tedious but straightforward, and the results are 



(7) 
(8) 



Mi = -((if - UbfPb), 

+ (p f - UbfPbf + Uifttbjbibtbi) - pi), (9) 



(4 



3t* 

2~(M/ - UbfPb) - (Pf - UbfPbf 



ipfU b f 

x {{b\hb\b t ) - pi) + Uifdbtbibtbiblbi) - pi), (10) 

where the operator averages (b\bib\bi- ■ •) can be easily 
calculated from the knowledge of w nb as ((b\bi) k ) — 



Pb 



w„ b n b 



(11) 



n b =0 



is the filling of the bosons. Therefore, the accuracy of the 
calculations can be further checked by comparing the real- 
axis calculation of the moments p^ by directly performing 
the integrations over the real frequency to the exact results 
given in Eqs. (0 - ( TTOb . 

Similarly, the spectral moment sum rules for 
the retarded self-energy are defined as = 
-(l/7r)/^ eiwa; m ImE(a;). One can also calculate these 
moments by using Dyson's equation, and after some 
significant algebra, the results are 



C$ = U b 2 f((b\b i b\b i )-pl), 
Of = -Ulf(p s + 2U b fPb)({b\hb\b l 
+ Ulf{{b\hb\b l b\b l )-pl). 



Pi) 



(12) 



(13) 



In this way, the accuracy of the calculations can be again 
benchmarked by comparing the real-axis integration of the 
spectral moments with the exact values given in Eqs. (fT~2t 
and ( fT3l l. In addition, we would like to mention that the large 
frequency limit of the real-axis self-energy 



E(w — > 00) = UbfPb 



(14) 



provides another independent check of the numerics. 

Typically, there are a number of poles in the self-energy 
near the correlation induced band edges when Ubf is large 
enough, and in order to guarantee the accuracy of the sum 
rules that are calculated on the real-axis by direct quadra- 
ture, the pole contributions have to be included by hand as 
described here. Using Eqs. <J2J and (0 with iu> n — > uj + i0 + , 
we find that the locations ui p of these poles are determined by 
the transcendental equation 



E 

n b =0 



'Jn b YL ( uj p + ^ U bf m b ) 



= 0. 



(15) 
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In a typical calculation, we restrict the bosonic occupancies 
to be less than some maximal number of multiple boson oc- 
cupancy. Then the transcendental equation becomes a fi- 
nite equation which is easy to solve directly. Using such a 
procedure, the self-energy can then be written as = 
S rcg (w) + J2 P Rpl (w - w P + «0 + ), where £ rcg (w) is the reg- 
ular piece, p sums over the finite number of poles that were 
found to solve the truncated equation, and R p is the residue 
of the pole ui p that is calculated next (only poles with positive 
residue are included in the sum as described below). Since 
the pole contribution dominates the spectral functions near the 
poles, one can expand Eq. (0 with ioj n — > ui + i0 + in powers 
of u> + fXf — and obtain G(lu) « gi/[uj + [if — £(u;)] + 

gs/[uj + pf — S(oj)] 3 + • ■ • where g\ = f_ dep(e) = 1 and 
03 = dep(e)e 2 = t* 2 /2. Plugging this expansion into 
Eq. (0 to obtain near the poles, and after some algebra, 
leads to 



Rn 



*0nJIlm 6 =4n»( w P + Pf ~ U bf m b ) 



E^=0 W n b Y,n> b *n b Um b ^{n b ,n' b }( U P + M/ ~ U bf m b ) 



(16) 

where <5y is the Kronecker delta function. Only poles where 
the residue R p is positive correspond to real physical poles. 
As a result, the pole contribution to the self-energy moment 
Cm i s J2 P > Rp'^p, where p' sums over the poles with positive 
residues. 

We recall that the spectral moment sum rules given in 
Eqs. (0 - ( [Tol l and Eqs. ( TT2| > and ( fT3l ) have exactly the same 
form in the case of the Fermi-Fermi FK model fm l20ll2ill. 
However, in that case {{b\bi) k ) = (b\bi) = wi for any 
k due to the fermionic anticommutators. In addition, since 
w nb >2 — 0, the pole equation given in Eq. (fT5l l reduces con- 
siderably leading to a single pole at lj% = —p,f + U b f(l — Wx), 
for which the residue equation given in Eq. ( fT6] i reduces to the 
well-known result Ri = u>i(l — wi)U b j — g% 112711 . 



C. Asymptotic Behavior for Large Frequencies 

There is another important application of the spectral mo- 
ment sum rules given above. They can be used to evaluate the 
high-frequency asymptotic behavior of the Green's function, 
self-energy and dynamical mean-field exactly. Therefore, this 
knowledge can be used to reduce the number of Matsubara 
frequencies that is needed to solve the imaginary-axis equa- 
tions ||2lll . 

When the Matsubara frequency is high enough, it can be 
shown that the Green's function can be written as G(iu> n ) — 
J2m=o A*m/(* w ™) m+1 ' an d similarly the self-energy can be 
written as E(tw„) = E(oo) + Erf C ™/K) m+ ' where 
£(oo) is a real constant, and it is the large-frequency limit 
of the self-energy. These expansions follow from the defini- 
tion of the spectral moment sum rules and the spectral for- 
mula for the retarded Green's function and self-energy. Plug- 
ging these expansions into Eq. (0, and using the definition 



of G 1 (iLu n ) given below Eq. (0, we obtain the asymp- 
totic expansion of the dynamical mean-field as X(iu> n ) — 
t* 2 /[2{iu) n )]-t* 2 (n f - U bf p b )/[2{iu n ) 2 ] + ■ • • . This asymp- 
totic expansion along with the expansion above allow us to 
treat the high-frequency tails of some quantities as described 
next. 

Let's take an energy cutoff e c which is much larger than the 
bandwidth of the interacting DOS, and use the asymptotic ex- 
pansion to sum over the Matsubara frequencies that are higher 
than e c . This also defines a cutoff n c for the Matsubara fre- 
quencies, given by the closest one to e c but lying below it, i.e. 
u) nc — (2n c + 1)tt//3 < e c . Plugging the asymptotic expan- 
sions into Eq. (0, we obtain 



2 exp I (3 



M/ ~ U bf n b U bb 

h fi b n h —n b (n b - 1) 



n 

n=0 



G Q x {iu) n ) - U bf n b 



iu> 

i*2 



t* 2 {nf-U bfPb ) 



n 

n—n c -\-l 
2 



1 M/ - Ub/n b 



2(iw„) 3 



(17) 



Here, since the complex conjugate of iio n is equal to £u;_ n _i, 
we take the negative Matsubara frequencies into account by 
writing the absolute value squares in the products. To ap- 
proximate the infinite products that have an infinite number 
of terms, we first rewrite the infinite products as the exponen- 
tial of the sum of the logarithm of the individual terms. Then, 
for temperatures much lower than the bandwidth of the inter- 
acting DOS, we replace the sum by an integral, and then con- 
vert the integral over frequency to an integral over z = l/u>, 
leading to 



n 

n=n c -\-l 



i 



iuJ n (iuj n ) 2 {iuJ n ) 3 



exp ■ 



2tt 



I % In [(1 + bz 2 ) 2 + z 2 (a - cz 2 ) 2 } }. (18) 
'o 2 ' 



Here, a = fif — U b fn b , b — —t* 2 /2, and 



U b fp b )/2. As a result, the use of asymptotic expressions for 
the Green's function, self-energy, and the dynamical mean- 
field allow us to reduce the computational effort considerably 
for the imaginary-axis calculation by keeping e c small. Be- 
cause the DOS extends farther out in energy as T rises, one 
needs to increase e c in order to achieve the same level of ac- 
curacy for high T. 

The asymptotic expansions can also be used to treat the 
tails in the summation for the fermion filling, i.e. pf — 
0-/P) J2n=-oo G ( iuj n), which leads to 



PS = — 



48 
Mi 



1 n c — l 



(iuj n ) 2 (iw„) 3 (iuJ n ) 4 



(19) 



Here, we use the standard representation of the Fermi-Dirac 
distribution function F(x) = (1//3) J2^L-oo l/(«w„ — x) and 
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its derivatives with respect to x when x — * 0, to evaluate the 
coefficients of the Y^=~oo l/(* a7 «) fe tv P e - Since the 

contribution from the summation (which needs to be evaluated 
numerically) vanishes rapidly for large cj n , the filling can be 
calculated quite accurately. 



IV. NUMERICAL RESULTS 

Having described the DMFT formalism for the Fermi-Bose 
FK model, next we present some of our numerical results that 
are obtained by solving Eqs. ©, (0 and © self-consistently 
for given fermion and boson fillings, pj and pb, respectively. 
For this purpose, we choose a large energy cutoff e c = 2004*, 
and treat the tails of the high-frequency summations using the 
formalism described above. Since the system is unstable for 
Ubb < 0, and the Ubf < case can be mapped onto the 
Ubf > case as discussed in Sec. Ill Bl we restrict our anal- 
ysis to {Ubb, Ubf} > values. We remark here that numeri- 
cal solutions with w nb ^oo 7^ are unphysical, indicating an 
instability in the system corresponding to a bosonic collapse 
when Ubb < 0. 



A. Mapping onto the Spinless Fermi-Fermi FK Model 



We argued in Sec. Ill B I that the Fermi-Bose FK model can 
be exactly mapped onto the well-studied spinless Fermi-Fermi 
FK model for all parameter space at T = as long as the mix- 
ture is thermodynamically stable. Here, we numerically show 
how the system evolves away from the Fermi-Fermi behavior 
at higher temperatures. 

In Fig.Q] we calculate the region where the total probability 
of finding a boson in just two of the rib states is higher than 
99%. Therefore, below each line in all figures, the system 
can be effectively mapped onto the spinless Fermi-Fermi FK 
model. In Fig.QJa), we show temperature T vs. boson-boson 
repulsion Ubb mapping diagram for different boson fillings pb 
when bosons and fermions are uncoupled, i.e. Ubf = 0. As 
Pb increases, it is seen that the mapping is possible only for 
lower T values at a given Ubb, and is possible only for lower 
Ubb values at a given T. This is because it is energetically 
more favorable to have an occupation of multiple rib states as 
a function of increasing T and/or decreasing Ubb due to the 
Bose distribution function. 

When Ubf ^ 0, we show a T vs. Ubb mapping diagram for 
different Ubf in Fig. [Tib), and a Ubf vs. Ubb mapping diagram 
for different Tin Fig. [Tfc). As Ubf increases, it is seen that the 
mapping is possible for smaller and smaller parameter space 
compared to the Ubf = limit. This is because the coupling 
between bosons and fermions induces an attractive interaction 
between bosons such that the effective boson-boson repulsion 
Ufo decreases by some amount that is proportional to U^j to 
the lowest order in Ubf - Therefore, increasing Ubf leads to an 
occupation of additional rib states as discussed next. 



B. Density of States (DOS) for the Fermions 

In this subsection, we present our numerical results for the 
probability w nb of finding a boson in state rib, and the single- 
particle many-body density of states (DOS) for the fermions, 
to illustrate typical properties of the Fermi-Bose FK model. 
The probabilities are shown in Fig. [2 a) as a function of T for 
the first five boson occupancies when pb = 1, pf = 0.25, 
Ubb = t*, and Ubf — l.7t*. It is seen that w nb — 5\ nb at T = 
which is due to p b = J2n b =o w n b nb with £,T 6 =o w «t, = 
However, the occupation of the rib = 1 state decreases at finite 
T, while that of higher and higher rib states become finite with 
increasing T. 

The occupancy of multiple rib states has a strong effect 
on the many-body DOS for the fermions, i.e. A(w) = 
— (l/ir)ImG(iu> n — ► oj + ?'0 + ), and is given by 



A{w) 



1 



-Im 



dep(e) 



+ // / -E(w)-e + i0+ ! 



(20) 



where p(e) is the noninteracting DOS for the infinite- 
dimensional hypercubic lattice defined below Eq. (O, and the 
infinitesimal is needed only when ImS(cj) = 0. In Fig. [2b), 
there is a single Gaussian peak at T = corresponding to 
rib = 1 (which actually becomes a noninteracting system), but 
there are three peaks at T — O.lt* corresponding to rib = 0,1 
and 2. In addition, two more peaks occur at T = 4* cor- 
responding to rib = 3 and 4. This is again because it is 
energetically more favorable to have an occupation of mul- 
tiple rib states as a function of increasing T due to the bosonic 
character of the heavy atoms. We mention that while there 
is not any pole in the self-energy when T = 0, there are 
two physical poles at u>\ w — 0.348i* and u>2 ~ 2.1494* 
with residues R x w R 2 « 0.166i* when T = O.lt*, and at 
lji w 0.7254* and u 2 ~ 2.8714* with residues R t sa 0.4194* 
and i? 2 w 0.4024* when T = 4*. Notice that poles occur near 
the correlation induced band edges. 

We investigate the effects of Ubf on A(uj) in Fig. [2jc), 
where it is seen that increasing (decreasing) Ubf increases 
(decreases) the number of peaks. This is again because the 
coupling between the bosons and fermions induces an attrac- 
tive interaction between bosons such that the effective boson- 
boson repulsion decreases, which leads to an occupation of 
additional rib states. For instance, as shown in the figure, there 
is a single peak when Ubf = 4*, but there are three of them 
when Ubf = 1.74*. In addition, the band gap between differ- 
ent rib states increases with increasing Ubf. We again men- 
tion that while there is not any pole in the self-energy when 
Ubf = 4*, there are two physical poles at u>i w —0.3484* 
and w 2 « 2.1494* with residues Ri w R 2 « 0.1664* when 
U b f = 1.74*, and at wi w 0.06054* and uj 2 ~ 4.2954* with 
residues Ri w R 2 w 1.7594* when Ubf = 34*. Notice again 
that poles occur near the correlation induced band edges. 

Although the DOS for the fermions shows rich structures, 
this is not yet a measurable quantity in ultracold atomic sys- 
tems (it might be feasible with an appropriately designed rf- 
frequency experiment that acts like a photoemission experi- 
ment of the fermions 12811 . but there one would be observ- 
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FIG. 1: Mapping diagrams show when the Fermi-Bose FK model can be mapped onto the spinless Fermi-Fermi FK model as a function of 
boson-boson repulsion Ubt- In <*H figures, mapping boundaries separate the regions where total probability of finding a boson in two of the ni 
states is higher than 99% (which occurs below each line). 




FIG. 2: (a) The probability w nb of finding a boson in state rib is shown as a function of temperature T for the first five boson occupancies, (b, 
c) The single-particle many-body density of states for the fermions A(uj) is shown as a function of frequency uj. 



ing the DOS multiplied by the Fermi-Dirac distribution func- 
tion, the so-called lesser spectral function). For this reason, 
we next discuss the momentum distribution of the fermions 
and bosons which could be easily measured in a time-of-fiight 
measurement. 



C. Momentum Distribution 

The occupancy of multiple rib states has also significant 
effect on the momentum distribution of the fermions, i.e. 
n m d(k) = (l//?)I]^L_ 00 G(k, iw„), which becomes 

1 °° 1 
rw(k) = - V — — jt-, (21) 

where G(k, iu> n ) is the momentum-resolved Green's function 
and e(k) = —2tfJ2i=i cos(kia) is the energy dispersion for 



the fermions defined below Eq. ©. The momentum distribu- 
tion of the bosons is simply a constant given by the filling of 
the bosons since they are localized with zero tunneling. To 
regularize the frequency summation in our numerical calcula- 
tions, we subtract (1/(3) X^=-oo l/[^n+M/- s (°°)-e(k)], 
and add F[—fif + S(oo) + e(k)] in Eq. (ED, where F(x) is 
the Fermi-Dirac distribution function. 

Our numerical results for n m( j(k) vs. the fermion disper- 
sion e(k) are shown in Fig. [3] As shown in Fig. [3ja), n me j(k) 
is a step function at T = 0, and it broadens as a function 
of increasing T, which is expected due to Fermi-Dirac statis- 
tics. For a fixed T, the effects of Ubf on n m d(k) are shown 
in Fig.|3jb), where it is seen that increasing Ubf also broad- 
ens n mc i(k) just like the temperature, which is the expected 
many -body effect due to a finite lifetime of the fermionic exci- 
tations. In the FK model, the fermions are usually not a Fermi 
liquid (or more correctly a Fermi gas in our context) at T = 0, 
but for the case of integer boson fillings pb = 0, 1, 2, • • • , oo 
(presented here) they are, because the system evolves to an 
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o 1 ■ ~~ 

-16 -8 8 16 

e(k)/t* 



(b) p b - 1, pf - 0.25, U bb = t*, and T = 0.1t* 




-16 -8 8 16 

e(k)/t* 

FIG. 3: The momentum distribution of the fermions n m d(k) is 
shown as a function of the dispersion e(k). The parameters used 
in these figures are the same as the ones used in Fig. [2] 



effective noninteracting system as T — > 0. 



D. Average Kinetic Energy 

Another important quantity that can be measured in ultra- 
cold atomic systems is the average kinetic energy of the par- 
ticles, which is given by (e(k)) = e(k)n m£ ;(k) for the 
fermions, and vanishes for the bosons since they are local- 
ized with zero tunneling. Since the momentum distribution is 
measured directly in experiment, one can process the data to 
determine the average kinetic energy in the lattice under the 
assumption that the momentum distribution has not changed 
significantly during the time-of-flight experiment. For our cal- 
culations, we convert the summation over momentum to an 
integral over energy. The resulting expression is 

/oo 
dep(e)en md (e), (22) 
-oo 

where p(e) is the noninteracting DOS for the infinite- 
dimensional hypercubic lattice defined below Eq. ©. 

Our numerical results for (e(k)) vs. T are shown in Fig. [4] 



-0.24 1 1 1 1 1 

0.5 1 1.5 2 

T/t* 

FIG. 4: The average kinetic energy of the fermions (e(k)) is shown 
as a function of temperature T. 



When Ubf is small, e.g. Ubf = t*, it is seen that (e(k)) in- 
creases monotonically as a function of T. However, when 
Ubf becomes large enough, see e.g. Ubf = 1.7t*, it has a 
local minimum at finite T after an initial increase for lower 
temperatures. This minimum moves towards higher T as Ubf 
becomes larger, compare e.g. Ubf — 1.7t* case with that of 
Ubf = 3t*. The nonmonotonic behavior occurs only for large 
enough Ub / values and is a consequence of the correlation in- 
duced band gaps that are present in the many-body DOS, as 
discussed in Sec. IIV Bl This may be one of the most direct 
ways, with current available technology, to infer the changes 
in the DOS as a function of T in experiment. 



V. CONCLUSIONS 

In this work, we analyzed Fermi-Bose mixtures consisting 
of light fermions and heavy bosons that are loaded into op- 
tical lattices. To describe such mixtures, we considered the 
Fermi-Bose version of the FK model, the Fermi-Fermi version 
of which has been widely discussed in the condensed matter 
literature. In our model, we assumed that the bosons are local- 
ized such that their tunneling to other sites vanishes (i& = 0) 
but that the system can statistically sample all low energy con- 
figurations of the heavy atoms. This perspective makes sense 
in ultracold atomic experiments if the difference in the tunnel- 
ing amplitudes between the light fermions and heavy bosons 
is large enough so that the quantum nature of the bosons can 
be neglected, but they can reorganize their positions to allow 
the system to sample different configurations. An alternative 
perspective is that once the optical lattice has been introduced, 
the heavy bosons become frozen into a specific configuration, 
that is randomly chosen from the configurations that are ener- 
getically favorable in the statistical mechanical ensemble; as 
one repeats many experiments and averages over the different 
configurations, one would then reproduce the results of the FK 
model described here J29ll . 

First, we discussed the symmetries of the Hamiltonian, and 
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showed that the Fermi-Bose FK model can be mapped exactly 
onto the spinless Fermi-Fermi FK model at zero temperature 
for all parameter space as long as the mixture is thermody- 
namically stable. Since this mapping is only approximate at 
low temperatures and fails at high temperatures, we developed 
DMFT to investigate the effects of temperature (recall that the 
DMFT becomes exact in infinite dimensions). We also cal- 
culated spectral moment sum rules for the retarded Green's 
function and retarded self-energy, and used them to check the 
accuracy of our numerical calculations, as well as to reduce 
the computational cost. When bosons and fermions are un- 
coupled (Ubf — 0), we showed that, as the boson filling in- 
creases, the mapping is possible only for lower T values at a 
given Ubb, and is possible only for lower Ubb values at a given 
T. This is because it is energetically more favorable to have 
occupation of multiple nj, states as a function of increasing T 
and/or decreasing Ubb due to the Bose statistics. As Ub / in- 
creases, we found that the mapping is possible for smaller and 
smaller parameter space compared to the Ubf = limit. This 
is because the coupling between bosons and fermions induces 
an attractive interaction between bosons such that the effec- 
tive boson-boson repulsion decreases by some amount 
that is proportional to to the lowest order in Ubf. There- 
fore, increasing Ubf leads to an occupation of additional rib 
states. 

We also presented typical numerical results for the Fermi- 
Bose FK model including the occupancy of bosonic states, 
single-particle many-body DOS for the fermions, experimen- 
tally relevant momentum distribution, and the average kinetic 
energy. We found that the occupancy of multiple bosonic 
states has a strong effect on the DOS for the fermions, leading 
to strong modulations as a function of frequency. The number 
of peaks corresponds to the number of bosonic states that are 



occupied, and it increases as a function of increasing T and/or 
Ubf. In addition, we showed that increasing Ubf at a fixed 
T broadens the momentum distribution of the fermions, just 
like the effects of temperature by itself. We also showed how 
the average kinetic energy evolves with T and how one can 
infer changes to the DOS via structure in the average kinetic 
energy. 

We hope that some of these results could be experimentally 
realized in ultracold atomic systems. We think, for instance, 
K-Rb, Li-K, or Li-Cs mixtures are a good initial candidates 
for simulating the Bose-Fermi FK model. In addition, one 
could also create species-dependent optical lattices for differ- 
ent isotopes of the same atom such that the bosonic isotope is 
localized but not the fermionic one. 

Finally, if one recalls the local-density approximation as a 
first approximation to the effects of the trap in a real exper- 
imental system, we expect that the methods described here 
could be quickly used to generate approximate results for den- 
sity distributions across the trap and for a variation of the den- 
sity of states or of the momentum distribution. Such results 
are beyond the scope of this work, but could be investigated 
if one is interested in directly modeling a specific experiment 
that is described by the Fermi-Bose FK model. 
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